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Abstract 

We  consider  a diffuse-interface  model  for  fluid-fluid  systems.  In  classical  models, 
an  interface  between  two  fluids  is  treated  as  infinitely  thin,  or  sharp,  and  is  endowed 
with  properties  such  as  surface  tension.  Diffuse-interface  theories  replace  this  sharp 
interface  with  continuous  variations  of  an  order  parameter  such  as  density  in  a way 
consistent  with  microscopic  theories  of  the  interface.  Surface  tension  effects,  for  exam- 
ple, are  incorporated  into  the  model  through  a modified  stress  tensor  in  the  classical 
Navier-Stokes  equations.  We  relate  the  diffuse-interface  model  to  classical,  sharp  in- 
terface models  by  deriving  asymptotically  the  governing  equations  and  the  associated 
boundary  conditions  used  in  the  sharp-interface  formulation.  We  illustrate  the  diffuse- 
interface  approach  by  modehng  internal  gravity  waves,  which  have  been  observed  ex- 
perimentally by  Berg  et  al.  in  xenon  near  its  critical  point.  We  obtain  static  density 
profiles,  compute  internal  wave  frequencies  and  compare  with  their  experimental  data 
and  theoretical  (classical)  results  both  above  and  below  the  critical  temperature.  The 
results  reveal  a singularity  in  the  diffuse-interface  model  in  the  limit  of  incompressible 
perturbations. 
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1.  Introduction 


Classical  models  of  fluid-fluid  systems  treat  the  interfacial  region  between  two  fluids  as 
an  infinitely  thin,  or  sharp,  dividing  surface  and  endow  it  with  properties  such  as  surface 
tension.  Properties  such  as  density  take  on  their  bulk  values  all  the  way  up  to  the  point 
of  discontinuity.  The  equations  of  motion  are  solved  in  separate  domains  and  the  interface, 
where  appropriate  boundary  conditions  are  applied,  must  be  tracked  explicitly. 

In  diffuse-interface  theory,  the  interfacial  region  is  represented  by  continuous  variations 
of  an  order  parameter  such  as  density  in  a way  consistent  with  microscopic  theories  of  the 
interface  (e.g.  Rowlinson  and  Widom^).  These  density  variations,  which  give  the  interface 
nonzero  thickness  and  internal  structure,  connect  smoothly  to  the  bulk  values  of  the  density 
on  either  side  of  the  interface.  Interfacial  excess  quantities  are  defined  in  terms  of  the 
density  variations  through  this  region.  Consequently,  these  excess  quantities,  rather  than 
being  defined  on  a two-dimensional  surface,  are  distributed  throughout  a three-dimensional 
layer.  The  equations  of  motion,  modified  to  account  for  the  presence  of  this  layer,  apply  over 
the  entire  domain.  The  interfacial  region  is  then  identified  by  a range  of  constant  density 
contours;  no  interface  tracking  is  required. 

Diffuse-interface  theories  have  been  used  successfully  in  modeling  solidification  and  phase 
transitions.  A notable  example  is  the  computation  of  dendritic  growth  (e.g.  Kobayashi^, 
Wheeler  et  al.^  and  Warren  and  Boettinger^).  These  models  have  focused  on  diffusive  mech- 
anisms (thermal  and  solutal)  and  have  not  addressed  advective  transport. 

Recently,  there  have  been  studies  of  a variety  of  complex  hydrodynamic  phenomena  using 
the  diffuse-interface  approach.  Antanovskii®  developed  a diffuse-interface  model  describing 
fluid  flow,  with  heat  and  species  transport  for  a binary  fluid  where  composition  was  the  or- 
der parameter.  A capillary  tensor,  which  accounted  for  capillary  forces  associated  with  the 
interfacial  region,  was  derived  using  reversible  thermodynamic  arguments  and  the  surface 
tension  given  in  terms  of  the  excess  internal  energy.  This  model  was  illustrated  by  computing 
thermocapillary-driven  flow  in  a gap.  A similar  diffuse-interface  model  describing  the  motion 
of  an  isothermal  binary  fluid  was  considered  by  Gurtin  et  al.^,  whose  derivation  was  based  on 
microforce  balance  laws.  The  model  consisted  of  modified  Navier-Stokes  equations,  which 
included  a capillary  tensor  accounting  for  interfacial  forces,  and  a coupled  Cahn-Hilliard 
equation  for  the  composition  modified  to  account  for  hydrodynamic  transport.  A model 
for  incompressible,  immiscible  fluids,  where  the  order  parameter  was  a conserved  quantity 
transported  with  the  fluid  (e.g.  density),  was  also  described.  Those  authors  illustrated  their 
approach  by  computing  the  coarsening  and  the  associated  fluid  motion  for  a binary  fluid. 
Chella  and  Vinals^  used  this  model  to  compute  the  mixing  and  interfacial  stretching  of  a 
binary  fluid  in  a shear  flow.  Jasnow  and  Vinals®  studied  thermocapillary  drop  migration  and 
coalescence  as  well  as  spinodal  decomposition  using  this  model.  In  addition,  they  presented 
a derivation  of  the  model  using  a Hamiltonian  formalism.  Jacqmin^  addressed  a number  of 
complex  flows,  including  droplet  breakup,  wave-breaking  and  contact-line  motion  using  the 
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diffuse-interface  approach.  His  formulation  included  a wall  potential  to  model  the  surface 
energy  of  the  container  boundary.  All  of  the  above  computations  demonstrate  the  unique 
capabilities  of  the  diffuse-interface  approach  to  model  flows  involving  complex  interface  mor- 
phologies and  topological  changes. 

A situation  in  which  diffuse-interface  theory  arises  naturally  is  in  the  description  of  a 
fluid  near  its  critical  point.  Even  in  the  absence  of  complex  interface  shapes  as  described 
above,  this  situation  can  involve  nontrivial  interfacial  phenomena.  At  temperatures  below 
the  critical  point,  two  distinct  fluid  phases,  separated  by  an  interface,  coexist.  At  temper- 
atures above  the  critical  point,  only  a single  fluid  phase  is  present;  the  interface  ceases  to 
exist.  As  the  critical  temperature  is  approached  from  below,  the  interface  becomes  infinitely 
diffuse.  There  has  been  a wealth  of  literature  devoted  to  the  study  of  the  dynamics  of  a 
diffuse  interface  in  the  context  of  critical  phenomena.  FelderhoH°  derived  a set  of  equations 
governing  the  dynamics  of  the  diffuse  interface  near  the  critical  point  of  a pure  fluid  using 
a Lagrangian  formalism.  Langer  and  Turski^^  developed  a similar  hydrodynamic  model  to 
describe  condensation  (nucleation)  of  a vapor  near  its  critical  point  using  a coarse-grained 
(diffuse-interface)  approach.  Extensive  analysis  using  renormalization-group  techniques  have 
been  performed  on  a diffuse-interface  model  (‘Model  H’,  as  it  is  commonly  known  in  the 
literature^^)  which  describes  the  dynamics  of  a binary  fluid  phase  transition  as  well  as  a 
pure  fluid  near  its  critical  point  (e.g.  Halperin  et  al}^,  Siggia  et  and  Hohenberg  and 

Halperin^^).  A similar  modeH^  has  been  used  to  study  the  dynamics  of  a near-critical  fluid 
in  a shear  flow  (Onuki  and  Kawasaki^®  and  Onuki  et  al}^) 

These  diffuse-interface  models  share  common  features  with  those  developed  from  a more 
computational  point  of  view.  Brackbill  et  al}^  developed  a ‘continuum  surface  force’  model 
wherein  they  identified  a volume  force  which  represents  surface  tension  spread  over  a small 
but  finite  three-dimensional  interfacial  domain.  This  volume  force  was  related  to  a ‘color’ 
function  which,  for  example,  can  represent  density  for  incompressible  flows.  The  defining 
characteristics  of  this  volume  force  were  that  it  gives  the  correct  surface  force  in  the  limit 
of  a sharp  interface  and  is  nonzero  only  in  the  interfacial  region.  Another  related  approach 
has  been  described  by  Unverdi  and  Tryggvason^®’^°.  Their  approach  is  a front-tracking 
technique  which  employs  a numerically-diffuse  description  of  the  interface.  They  construct 
an  indicator  function,  based  on  the  known  position  of  the  (sharp)  interface,  which  identifies 
fluid  properties  such  as  density  and  viscosity.  This  function  is  then  artificially  spread  out 
over  a small  region  on  the  scale  of  the  computational  mesh  size,  allowing  the  fluid  properties 
to  vary  smoothly  through  this  interfacial  region.  The  surface  force  (i.e.  surface  tension)  is 
also  distributed  over  this  interfacial  region  so  that  a single-domain  approach  can  then  be 
used  to  calculate  the  flow.  This  flow  then  determines  how  the  interface  is  advected. 

Recently,  a unique  hydrodynamic  phenomena  that  occurs  near  a fluid’s  critical  point  has 
been  described  and  analyzed  by  Berg  et  They  found  that,  owing  to  its  large  compress- 
ibility near  the  critical  point,  xenon  stratifies  under  its  own  weight  and  that  this  density 
stratification  supports  internal  gravity  waves  when  perturbed.  In  contrast  to  geophysical 
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applications  where  length  scales  can  be  on  the  order  of  miles,  the  stratification  in  near- 
critical  xenon  occurs  on  the  scale  of  centimeters.  Experimentally,  they  measured  internal 
gravity  wave  frequencies  as  a function  of  temperature,  both  above  and  below  the  critical 
temperature.  Their  theoretical  development  consisted  of  two  separate  classical  models;  one 
that  applied  above  the  critical  temperature  (single-phase  region)  and  another  that  applied 
below  the  critical  temperature  (two-phase  region).  In  each  case,  they  used  the  restricted 
cubic  modeP^’^^  as  the  equation  of  state. 

It  is  in  the  context  of  critical  phenomena  that  we  shall  apply  the  diffuse-interface  approach 
to  fluid  motion.  Our  first  objective  is  to  formulate  a diffuse-interface  model  for  a pure  fluid 
near  its  critical  point.  Of  specific  interest  here  will  be  the  interpretation  in  terms  of  the 
underlying  global  quantities  (mass,  energy  and  entropy)  and  their  associated  global  balance 
laws.  Our  second  objective  is  to  analyze  the  model  in  the  sharp-interface  limit  in  order  to 
recover  both  the  governing  equations  and  interfacial  conditions  of  the  classical  model.  Our 
third  objective  is  to  illustrate  the  diffuse-interface  approach  using  a model  for  internal  gravity 
waves  in  xenon,  which  follows  the  work  of  Berg  et  The  present  approach,  which  applies 
both  above  and  below  the  critical  temperature,  shall  be  compared  with  their  experimental 
and  theoretical  results.  The  focus  here  is  to  use  this  application  to  highlight  the  diffuse- 
interface  model  rather  than  to  improve  on  the  good  agreement  between  the  experiments 
and  the  classical  theory  used  by  Berg  et  al?^  Consequently,  we  shall  use  a simpler  equation 
of  state,  which  allows  an  analytical  solution  for  the  density  profile,  rather  than  the  more 
accurate  but  complicated  restricted  cubic  model  used  by  Berg  et  al. 

In  Section  2 we  formulate  the  diffuse-interface  model  which  describes  the  compressible, 
adiabatic  motion  of  a single-component  fluid  near  its  critical  point.  In  Section  3 we  derive 
the  classical  sharp-interface  governing  equations  and  associated  boundary  conditions  from 
the  diffuse-interface  model.  In  Section  4 we  outline  the  configuration  used  by  Berg  et 
to  study  internal  gravity  waves  in  near-critical  xenon.  In  Section  5 we  obtain  static  density 
profiles  using  an  approximate  van  der  Waals  equation  of  state.  In  Section  6 we  compute 
internal  gravity  wave  frequencies  associated  with  these  profiles  and  compare  with  the  results 
of  Berg  et  al.  A discussion  of  the  results  is  given  in  Section  7,  and  a conclusion  is  given  in 
Section  8.  We  include  in  the  appendices  a derivation  of  the  diffuse-interface  equations  which 
take  into  account  viscous  and  thermal  dissipation. 

2.  Formulation 

The  hydrodynamic  equations  governing  inviscid,  compressible  flow  of  a single-component 
fluid  near  its  critical  point  (e.g.  Felderhoff°,  Langer  and  Turski^^,  Hohenberg  and  Halperin^^, 
Gurtin  et  al.^,  Jasnow  and  Vinals®  and  Jacqmin®)  are  given  by 

^ + V(pu)  = 0,  (la) 
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p 


— Vp  — pgz  + KV  ■ T 


(lb) 


(J  + („.V)u)  = 

ds 

— + u-\/s  ■=  0,  (Ic) 

where  p is  the  fluid  density,  u is  the  three-dimensional  fluid  velocity  with  components 
(u,  Ujtn),  p — p{s,p')  is  the  thermodynamic  pressure  specifled  by  an  appropriate  equation 
of  state,  s(T,p)  is  the  entropy  per  unit  mass,  T is  the  temperature,  g is  the  gravitational 
acceleration,  z is  the  unit  normal  in  the  vertical  direction,  K is  the  gradient  energy  coefficient 
which  for  simplicity  we  assume  to  be  constant,  and  T is  the  capillary  tensor  given  by 

T=  (^>VV  + llV;Op)/-V^®Vp,  (2) 

where  I is  the  identity  matrix  and  0 is  the  tensor  (outer)  product.  Note  that  T has  the 
property 


V-T  = ^V(VV).  (3) 

The  key  difference  between  these  equations  and  the  classical  equations  describing  a two- 
fluid  system  is  the  presence  of  a capillary  tensor  T which  models  capillary  forces  associated 
with  the  interface.  Korteweg^^  was  the  first  to  use  such  a term  to  describe  capillary  effects  of 
a diffuse  interface.  The  presence  of  this  term  allows  these  equations  to  apply  over  the  entire 
domain,  including  the  interfacial  region;  no  interfacial  conditions  are  required.  Boundary 
conditions  applied  at  container  walls,  for  example,  are  required.  These  shall  be  discussed 
more  thoroughly  as  needed  for  the  present  work.  The  intersection  of  the  diffuse  interface 
and  a container  boundary  is,  in  general,  a moving  contact  line  and  is  not  treated  here  but 
has  been  discussed  to  some  extent  by  Cahn^®  and  by  Jacqmin®. 

In  order  to  understand  the  nonclassical  capillary  terms  appearing  in  the  momentum 
equation  and  in  the  local  entropy  balance,  it  is  helpful  to  relate  equations  (1)  to  global 
quantities  and  balance  laws.  We  define  the  mass  M,  energy  E and  entropy  5 in  a material 
subvolume  Q,[t)  of  the  total  volume  V by 


The  total  mass  is  given  simply  by  the  integral  over  the  local  density.  The  total  energy  is 
composed  of  classical  and  nonclassical  contributions.  The  classical  terms  include  kinetic 
energy,  gravitational  potential  energy  and  internal  energy  e[s,p)  given  per  unit  mass.  In 
addition  to  these  classical  contributions  there  is  a nonclassical  contribution  in  the  form 
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of  a gradient  energy  associated  with  steep  variations  in  density.^®  This  nonclassical  term 
represents  an  energy  excess  associated  with  the  interfacial  region.  Consistent  with  a constant 
gradient  energy  coefficient  in  equation  (4b)  there  is  no  gradient  entropy  term  in  equation  (4c) 
(e.g.  see  appendix  B). 

The  governing  equations  (1)  are  consistent  with  the  balance  laws 


dM 

dt 

d^ 

dt 

dt 


0, 


/ ( 

J6n{t)  \ 


Dp 

—pu  ■ h + Kn  ■ T ■ u -\-  K P ' ^ 


dA, 


(5a) 

(5b) 

(5c) 


where  6^n(t)  is  the  boundary  of  the  subvolume  and  n is  its  outward  normal.  The  first  balance 
law  (5a)  simply  represents  conservation  of  mass.  The  second  balance  law  (5b)  states  that 
the  change  in  the  energy  of  the  subvolume  is  associated  with  the  rate  of  (reversible)  work 
done  by  pressure  and  capillary  forces  on  the  boundary  as  well  as  a nonclassical  flux  of  stored 
energy  in  the  interface.  Note  that  this  flux  term  is  associated  with  compression  of  the  flow 
in  the  interfacial  region.  Wang  et  al?'^  identified  a similar  nonclassical  entropy  flux  term  in 
their  phase-field  model  of  solidification  (see  their  equation  (6)).  Whereas  their  term  involved 
the  partial  derivative  of  the  order  parameter  with  respect  to  time,  our  term  involves  the 
total  derivative,  Dpj Dt  = + it  • Vp,  since  we  have  accounted  for  fluid  motion.  They 

identified  this  term  as  an  entropy  flux  associated  with  variations  in  the  phase-field  at  the 
boundary  of  the  subvolume.  The  third  balance  law  (5c)  states  that  the  change  in  total 
entropy  in  the  subvolume  Q,{t)  is  conserved  (i.e.  zero  entropy  production).  Note  that  since 
we  have  neglected  dissipation,  these  equations  describe  an  adiabatic  process. 

A full  derivation  of  the  diffuse-interface  equations  that  account  for  dissipative  effects  and 
for  variations  in  the  gradient  energy  and  entropy  coefficients  is  given  in  Appendix  B. 


3.  Sharp-Interface  Limit 

Owing  to  the  nonclassical  nature  of  diffuse-interface  models,  there  has  been  much  effort  de- 
voted to  relating  these  models  to  their  sharp-interface  counterparts.  The  idea  behind  the 
diffuse-interface  approach  is  that  the  solutions  of  the  diffuse-interface  equations  represent 
asymptotically  the  solutions  of  the  classical,  sharp-interface  equations.  Extensive  asymp- 
totic analyses  have  been  performed  on  phase-field  models  of  solidification  to  show  that  they 
recover  a number  of  sharp-interface  solidification  models  (e.g.  see  Caginalp^®),  and  further 
comparisons  between  sharp-  and  diffuse-interface  models  of  solidification  can  be  found  in 
Braun  et  al?^  Diffuse-interface  models  of  fluid-fluid  interfaces  have  received  less  attention  in 
this  regard.  Antanovskii®  derived  from  the  diffuse-interface  model,  the  classical  hydrostatic 
balance  for  the  case  of  a flat  interface  in  equilibrium  and  the  Laplace-Young  equation  for  the 
case  of  a spherical  interface  in  equilibrium.  Jasnow  and  Vinals®  derived  from  the  capillary 
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term  in  their  momentum  equation  the  appropriate  sharp-interface  tangential  and  normal 
forces  when  the  surface  tension  was  a slowly  varying  function  along  the  interface. 

In  this  section  we  derive  from  equations  (1)  the  associated  sharp-interface  equations  and 
boundary  conditions.  A complete  description  of  the  sharp-interface  boundary  conditions 
appropriate  for  two-phase  systems  can  be  found  in  Delhaye^®. 

We  begin  by  writing  K = e^K  = e^[gL/ pc)K*  in  (1)  where  e is  a small  parameter 
measuring  the  thickness  of  the  interface  and  K*  is  dimensionless  and  0(1).  It  is  clear  that 
if  we  formally  take  the  limit  e — > 0 in  the  bulk  regions  away  from  the  interfacial  layer  we 
recover  at  leading-order,  the  classical  equations  appropriate  for  the  bulk  phases, 


Dp 

Dt 

= -pV  • u. 

(6a) 

Du 

(6b) 

Dt 

= -Vp-pgz, 

Ds 

(6c) 

'm 

= 0. 

Next,  we  seek  the  sharp-interface  boundary  conditions.  Since  the  diffuse-interface  gov- 
erning equations  apply  over  the  entire  domain,  including  the  interfacial  region,  we  shall  use 
standard  pillbox  arguments  to  obtain  these  conditions.  We  define  the  cylindrical  pillbox  as 
follows.  We  consider  the  contour  of  density  upon  which  the  interfacial  region  collapses  in 
the  limit  e — > 0.  The  pillbox  (shown  in  figure  1)  encloses  a portion  of  this  surface  at  a fixed 
point  in  time  in  such  a way  that  the  top  of  the  pillbox  is  above  the  surface  at  a height  ( = S 
and  the  bottom  of  the  pillbox  is  below  the  surface  at  a height  (=  —6.  Here,  ( is  a local 
coordinate  which  measures  distance  from  the  surface  in  the  normal  direction.  We  define 
normal  vectors  htop,  ribot  ^■iid  m for  the  top,  bottom,  and  side  of  the  pillbox,  respectively. 
The  pillbox  position  is  then  fixed  and  the  interface  is  allowed  to  sweep  through  at  the  next 
increment  of  time.  The  key  limit  in  the  pillbox  argument  is  that  e <C  ^ -C  T where  L is 
an  0(1)  length  scale  associated  with  the  outer  flow.  In  this  limit,  the  volume  of  the  pillbox 
becomes  negligible  on  the  outer  scales  but  the  variations  in  the  density,  which  define  the 
interfacial  region,  occur  over  a region  fully  contained  within  the  pillbox. 

Mass  balance:  Consider  first  the  mass  balance  (la).  We  integrate  this  equation  over 
the  volume  of  the  pillbox  Vp  to  obtain 


In  the  limit  e <C  ^ <C  T we  have  the  properties  that 

dp 


JV-o  Ot  JSr, 


f - f {pu)ui  ■ ndS, 

JVr,  ot  JS^ 


'vp  dt 


(7) 


(8a) 

(8b) 
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where  Uj  = (0,  O^wj^x^y))  is  the  velocity  of  the  surface  described  above  and  Sp  denotes  the 
surface  of  the  pillbox.  Note  that  V ■ uj  = 0.  It  then  follows  from  equation  (7)  and  the 
divergence  theorem  that 

0 = f p{u  — ui)  ■ nsdS,  (9) 

JSp 

where  ns  is  the  outward  normal  to  the  pillbox.  If  we  further  break  up  the  surface  integral 
into  contributions  from  the  top,  bottom  and  side  of  the  pillbox  we  find  that 


0 


I ^/)  ’ '^topdAiop  I pi^'^ 

^ Atop  At)ot 

+ y J p(u  — uj)  • md(dl, 


(10) 


where  we  have  written  the  surface  integral  over  the  side  of  the  pillbox  in  terms  of  a line 
integral  on  the  surface  and  an  integral  in  the  normal  direction.  In  the  limit  e <C  6'  <C  7/  the 
contribution  from  the  side  is  negligible  since  the  integrand  is  bounded.  Further,  in  this  limit, 
T^top  — '^hot  — ^ and  -A^op  — .^6ot  — A so  that 

0 — f p{u  — uj)  ■ h\'^  dA,  (11) 

J A 

where  indicates  the  jump  (top  to  bottom)  across  the  interface.  Since  the  interfacial  area 
A is  arbitrary,  we  must  have 


p{u  — uj)  • = 0. 


(12) 


This  is  the  classical  mass  balance  at  an  interface.  Note  that  consistent  with  the  original 
mass  balance,  this  boundary  condition  allows  mass  flux  across  the  interface. 

Momentum  balance:  Consider  next  the  momentum  equation  (lb).  If  we  integrate  over 
the  pillbox  as  before,  we  find  that 


" - k{ 


pgz  + V • [pu  (g)  u — pu  (g)  u/  + pJ  — KT]  \dV, 


(13) 


where  we  have  used  equation  (8b)  and  the  fact  that 


Du  d{pu) 
dt 


Dt 


+ V • {pu  0 u) . 


(14) 


In  the  limit  e S L,  the  gravitational  term  in  (13)  is  bounded  and  does  not  contribute. 
We  apply  the  divergence  theorem  to  the  remaining  terms  to  obtain 


0 


= / {pu  u — pu  <Si  uj  + pi  — KT)  ■ nsdS 

JSp 

= / {pu{u  - uj)  ■ hs  A pfis  — KT  ■ ns}  dS. 

J Sr. 


(15) 
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We  break  up  the  above  surface  integral  into  pieces  evaluated  on  the  top,  bottom  and  sides 
of  the  pillbox.  Note  that  in  the  limit  e <C  ^ <C  i/,  local  to  the  interface  we  have 


Vp^-n 


vV 


d( 


2 ’ 


(16) 

(17) 


so 


that 


n. 


m. 


(18) 

(19) 


Since  the  top  and  bottom  integrals  are  evaluated  at  ^ and  ( = —S,  respectively,  and 
the  variation  in  density  occurs  on  a scale  measured  by  e,  we  find  that  in  the  limit  e S L 
the  nonclassical  terms  do  not  contribute  to  the  top  and  bottom  surface  integrals.  This  leaves 
us  with  the  expression 


0 = J (^pu[u  — uj)  ■ h\'^ ph\'^^  dA 

^ J {pu[u  — Ui)  ■ rh  + pm  — KT  ■ m}  d(dl. 


(20) 


We  again  argue  that  the  variation  in  the  classical  terms  associated  with  the  velocity  occurs 
on  an  outer,  0{L),  scale  so  that  these  terms  do  not  contribute  to  the  side  integral  in  the 
limit  S L.  It  then  follows  that 


where 


0 = y (^pu[u  — Uj)  ■ dA  — y 'yiridl, 

pS 

7 = / {— p + i^m  • T • m} 


(21) 


~ 6 


(22) 


Here  = e<^  and  we  have  applied  the  limit  S/e  — > oo.  We  can  identify  the  expression  for 
7 as  the  excess  Kramer’s  (Grand  canonical)  potential  (see  below).  We  next  use  the  surface 
divergence  theorem  (e.g.  Weatherburn^^) 


y (f)rndl  = J V s(f>dA  — J (V5  • n)  (f)ndA, 


(23) 
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which  holds  for  a scalar  (f)  to  further  reduce  equation  (21).  Noting  that  the  surface  area  A is 
arbitrary  gives 

pu{u  — ui)  ■ n|^  + pn|t  = Vs7  — (V5  • fi)  7n.  (24) 


This  is  the  classical  (inviscid)  momentum  jump  condition  across  a fluid-fluid  interface.  Note 
that  the  first  term  represents  an  additional  normal  force  which  is  nonzero  whenever  there  is 
a jump  in  density  across  the  interface  and  a nonzero  mass  flux.  It  remains  to  confirm  that 
the  expression  for  7 given  by  (22)  can  be  interpreted  as  the  surface  tension. 

We  appeal  to  the  following  thermodynamic  arguments  to  obtain  independently  an  ex- 
pression for  the  surface  tension.  Consider  a nonclassical  extension  of  the  free  energy  defined 
by  the  functional 

^ = Jv  {pKp,  T)  + \k\Vp?)  dV.  (25) 

An  equation  governing  the  equilibrium  density  profile  can  be  obtained  by  minimizing  this 
free  energy  subject  to  the  constraint  of  constant  mass.  Therefore,  we  compute  the  variation 

0 = 6{l^(^pf+'^K\Vp\^-Xpjdv'^,  (26) 

where  A is  a constant  Lagrange  multiplier.  Variations  in  p lead  to 

0 - / + p/p  - KV^P  - A,  (27) 


where  A = where  is  a chemical  potential.  This  has  a first 

integral  whose  form  local  to  the  surface  as  described  above  can  be  written 


const.  = —\p  pf  — -K 


(28) 


The  thermodynamic  result  that  p = p^ fp  shows  that  const.  = 

Next,  we  consider  the  Kramer’s  potential  energy  fi  = JF  — and  define  its  excess 

by 


where  ao  represents  a Gibbs  dividing  surface  and  the  superscripts  (±00)  refer  to  the  far- 
fleld  values  on  either  side  of  the  interface.  If  we  now  use  the  equilibrium  condition  (27)  to 
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find 

dC) 

2 Uc, 

I)’ 

1(21' 
2 [d(, 


+ pX-  ^(-“)/<-“)  _ - p<-~>)  Uc 

j + - p(”)/(~)  - 

+ p<“)|<fC 

j +p‘”>|<iC  (30) 


where  we  have  used  A = Note  that  p is  the  solution  of  equa- 

tion (27);  gravity  does  not  play  a role  in  determining  the  density  profile  local  to  the  interface. 
We  can  combine  the  two  integrals  in  equation  (30)  to  find  that  is  independent  of  the 
position  of  the  dividing  surface.  We  have, 


£{_p  + pM  + jr(,le  + lJ|) 


(31) 


where  we  have  taken  = 0 without  loss  of  generality.  This  expression  is  identical  to  the 
previous  expression  (22)  for  7,  which  confirms  that  the  correct  surface  tension  appears  in 
the  sharp-interface  jump  condition  (24)  derived  from  the  diffuse-interface  model.  By  using 
equations  (27)  and  (28)  this  expression  can  be  rewritten  in  the  more  familiar  form 


(32) 


We  note  here  that  this  expression  relates  the  gradient  energy  coefficient  K (or  e)  to  the 
surface  tension,  which  in  principle  can  be  measured  experimentally. 

Entropy  balance:  Consider  next  the  entropy  balance  (Ic).  We  can  rewrite  this  using 
the  mass  balance  (la)  giving 

0 = ^ + V ■ (psu).  (33) 

Based  on  the  same  arguments  that  were  used  to  derive  the  jump  condition  for  mass,  we  find 
that  the  associated  jump  condition  for  entropy  must  be 


ps[u  — uj)  ■ h\'^  = J = 0, 


(34) 
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where  J is  the  mass  flux  across  the  interface.  This  condition  implies  that  the  jump  in  entropy 
must  vanish  when  there  is  nonzero  mass  flux  across  the  interface. 

Finally,  if  we  impose  the  condition  of  incompressibility  V • it  = 0,  we  find  that  this  leads 
to  the  condition  u ■ n = 0.  Together  with  (12)  it  follows  that 

■ n u~  ■ n = Ui  ■ n.  (35) 

4.  A Model  for  Internal  Gravity  Waves 

We  illustrate  the  diffuse-interface  approach  by  modeling  internal  gravity  waves.  This  appli- 
cation is  based  on  results  of  Berg  et  al?^ , who  observed  internal  gravity  waves  in  near-critical 
xenon.  In  their  experiments,  a mesh  paddle  inside  a small  sample  of  xenon  was  driven  at  a 
fixed  frequency,  stopped,  and  then  residual  internal  wave  frequencies  were  measured.  This 
was  done  for  a range  of  temperatures,  both  above  and  below  the  critical  temperature.  For 
the  case  where  the  cell  was  oriented  horizontally  (so  that  the  paddle  rotates  about  a verti- 
cal axis,  see  figure  2),  which  is  the  case  we  shall  consider  here,  they  observed  two  distinct 
internal  gravity  wave  modes.  They  also  modeled  these  waves  using  two  separate  classical 
models,  one  that  applied  in  the  one-phase  region  (above  the  critical  temperature)  and  an- 
other that  applied  in  the  two-phase  region  (below  the  critical  temperature).  They  used  as 
their  equation  of  state  the  restricted  cubic  modeF^’^^,  which  gave  a very  accurate  description 
of  the  static  density  profiles  as  the  temperature  varied.  Their  theoretical  predictions  for  the 
internal  gravity  wave  frequencies  agreed  well  with  those  observed  experimentally. 

The  model  configuration  consists  of  a rectangular  cell  with  dimensions  0 < a:  < a^, 
0 ^ ^ s.nd  —L<z<L  where  gravity  is  in  the  — z direction  (see  figure  2).  A paddle 

is  shown  in  the  figure  for  reference  purposes  but  no  paddle  is  present  in  the  mathematical 
model.  The  physical  dimensions  used  by  Berg  et  al.  correspond  to  = 7.6mm,  ay  = 38mm 
and  L = 9.5mm.  These  values  will  be  used  exclusively  here.  The  results  of  Berg  et  al. 
indicate  that  while  this  geometry  is  a simplification  of  the  actual  internal  geometry  of  the 
cell  used  in  the  experiments,  it  is  a reasonable  approximation  in  terms  of  identifying  the 
internal  wave  modes  and  frequencies. 

We  shall  model  the  motion  of  the  fluid  in  this  cell  with  the  diffuse-interface  approach 
described  in  section  2.  Our  analysis  of  these  equations  shall  proceed  as  follows.  We  first  seek 
static  density  profiles  (vertical  stratification  and  no  flow)  using  an  approximate  van  der  Waals 
equation  of  state.  We  then  perturb  these  density  profiles  and  identify  the  associated  natural 
internal  wave  frequencies.  We  shall  compare  the  present  results  with  the  experimental  data 
and  the  theoretical  predictions  of  Berg  et  al.  Additional  comparisons  shall  be  made  by 
computing  the  internal  wave  frequencies  with  the  sharp-interface  model  of  Berg  et  al.  but 
with  the  equation  of  state  used  in  the  present  work. 
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5.  Basic-State  Density  Profiles 


The  steady  basic-state  solution,  denoted  by  subscript  ‘O’,  has  zero  flow,  is  isothermal  and 
horizontally  uniform,  and  satisfies 


dz 


~Pog  + 


d^Po 
dz^  ' 


(36) 


This  equation  determines  the  static  density  profile  when  an  equation  of  state  p = p{T^  p)  is 
given.  Note  that  since  the  system  is  isothermal,  it  is  convenient  to  treat  the  pressure  as  a 
function  of  T and  p rather  than  of  5(T,  p)  and  p.  It  will  be  useful  to  work  with  a free  energy 
so  we  employ  the  thermodynamic  relation 


P = 


(37) 


where  / is  the  Helmholtz  free  energy  per  unit  mass.  This  allows  us  to  integrate  equation  (36) 
giving 


K 


d^Po 

dz^ 


djpf) 

dp 


-f-  -f  Co, 


(38) 


where  Cq  is  a constant  of  integration.  An  equivalent  derivation  of  equation  (38)  involves 
minimizing  the  total  energy  E (equation  (4b))  subject  to  constant  total  entropy  S (equa- 
tion (4c))  and  constant  total  mass  M (equation  (4a))  over  the  total  volume  V for  the  case  of 
no  flow.  Note  that  in  order  to  focus  on  the  critical  point,  we  have  assumed  that  the  average 
density  is  equal  to  the  critical  density  pc  so  that  M takes  on  the  fixed  value  p^V. 

We  shall  consider  a classical  van  der  Waals  equation  of  state  (e.g.  Stanley^^,  Callen^^)  as 
a reasonable  starting  point  for  our  analysis.  Despite  its  shortcomings  in  terms  of  predicting 
the  detailed  behavior  such  as  critical  exponents,  the  van  der  Waals  equation  of  state  provides 
a simple  description  of  the  qualitative  behavior  of  a fluid  near  its  critical  point.  The  van  der 
Waals  equation  of  state  is  given  by 


9RTc  2\ 


(39) 


where  R is  the  universal  gas  constant,  pc  is  the  critical  density  and  Tc  is  the  critical  temper- 
ature. 

If  we  use  equation  (37)  to  write  the  van  der  Waals  equation  of  state  (39)  in  terms  of  a 
free  energy,  we  And  that  it  has  a double-well  structure  below  Tc  and  a single-well  structure 
above  Tc-  We  shall  assume  a free  energy  per  unit  volume  of  the  general  form 


Pf 


PcCl{T)  + PcC2{T) 


+ Bq 


±T^(P_^\ 
2 T,  \ P.  ) 


(40) 
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where  Ci(T')  is  a temperature- dependent  integration  constant  whose  form  is  not  determined 
by  the  van  der  Waals  equation  of  state  (39)  and  C2(T)  is  related  to  Ci(T).  The  determination 
of  Ci(r)  requires  an  additional  thermal  equation  of  state  (e.g.  see  Callen^^).  Also,  the 
parameters  Bq  and  are  treated  as  adjustable.  This  reduces  to  the  van  der  Waals  form, 
expanded  locally  near  the  critical  temperature  and  density,  when  Bq  = ^pcRTc  and  = 4. 
Equations  (38)  and  (40)  lead  to  the  dimensionless  equation  for  the  static  density  profile 

= gz  + a^fp (41) 


where  p = {p  - pc)lpc,  f = {T  - Tc)ITc,  z = zjL,  = K pH BqB^  and  g = pcgLjBo.  We 
have  taken  p = pc  at  z = 0.  An  approximate  analytical  solution  to  equation  (41)  can  be 
obtained  in  the  limit  e «C  1-  For  T < 0 (i.e.  two-phase  region)  and  z < 0,  we  use  the  method 
of  matched  asymptotic  expansions  to  obtain 


P = 


Pout 


— ay  — 


T tanh 


(42) 


where  pout  corresponds  to  the  root  of  the  equation  0 = gz+a^Tp+p^  that  has  p > ay—T.  For 
T > 0 (i.e.  one-phase  region)  we  can  express  the  solution  in  terms  of  the  regular  expansion 


P — po  + 


(a2r  + 3p2)4 


(43) 


where  po  is  the  real  root  of  0 = pz  + a^Tpo  + Pq.  Note  that  the  density  profiles  are  antisym- 
metric about  z = 0.  Using  equations  (40),  (37)  and  s = —Jt  we  can  compute  expressions 
for  the  basic-state  pressure  po{z)  and  entropy  So[z).  Note  that  the  entropy  50(2)  depends  on 
Ci(T).  Further  details  of  this  are  discussed  in  section  7;  for  now  we  note  that  specification 
of  5o(a)  is  not  necessary  for  the  development  that  follows. 

Static  density  profiles  are  shown  in  figure  3 for  temperatures  both  above  and  below  the 
critical  temperature.  The  parameter  values  used  to  calculate  these  profiles  are  €dl  = 10““*, 
a = 4.85  and  g = 1.631  x lO-*  where  = e^/p.  The  parameters  and  p were  chosen 
to  fit  as  closely  as  possible  the  density  profiles  shown  in  Berg  et  al.  far  above  the  critical 
temperature  (see  their  figure  1).  The  dashed  curve  corresponds  to  the  density  profile  at 
the  critical  temperature.  Above  the  critical  temperature,  there  is  a single  stratified  phase. 
Below  the  critical  temperature,  the  fluid  separates  into  two  stratified  phases.  Again,  we  point 
out  that  while  these  profiles  do  not  embody  the  quantitative  behavior  in  terms  of  critical 
exponents,  they  do  capture  the  qualitative  behavior  of  the  density  near  the  critical  point. 

It  is  interesting  to  note  that  there  is  a significant  amount  of  stratification  which  occurs 
over  a relatively  small  length  scale.  While  this  stratification  often  plagues  those  seeking  to 
make  precise  measurements  of  physical  quantities  of  near-critical  fluids  such  as  xenon,  it  can 
also  be  seen  as  a unique  feature  through  which  phenomena  that  normally  occur  on  much 
larger  length  scales,  such  as  those  common  in  oceanography  or  atmospheric  sciences,  can  be 
studied  in  the  laboratory. 
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6.  Internal  Wave  Frequencies 

We  follow  the  analysis  of  Berg  et  and  seek  neutrally  stable  wave  modes  by  introducing 
perturbations  to  the  basic-state  solution  as  follows 


p 

TT 

o 

II 

+ 

^(2) 

cos(gxx) 

cos{qyy)e"‘^\ 

(44a) 

u 

= 0 

+ 

u(z) 

sm{q^x) 

cos(5yy)e^"‘, 

(44b) 

V 

- 0 

+ 

v{z) 

cos{qxx) 

sin(5y2/)e'"‘, 

(44c) 

w 

- 0 

+ 

w[z) 

cos{qxx) 

cos(5yy)e*‘^*. 

(44d) 

p 

II 

O 

+ 

p{z) 

cos  (5x2:) 

cos(5yy)e^‘^‘, 

(44e) 

s 

= -30(2:) 

+ 

s{z)- 

cos  ( 5x0;) 

cos(5yy)e^‘^‘. 

(44f) 

where  = "kj  ja^  and  qy  = -Kk/ay  are  the  wavenumbers  in  the  two  horizontal  directions  with 
integer  values  for  j and  k and  uj  is  the  frequency.  Note  that  we  are  interested  in  natural  wave 
modes  confined  to  the  box  and  therefore  have  a discrete  rather  than  continuous  set  of  wave 
vectors.  In  fact,  the  modes  we  describe  below  have  j = k = \.  The  form  for  the  velocity 
components  has  been  chosen  to  allow  u • n = 0 on  each  wall. 

As  a simplification,  we  shall  consider  incompressible  perturbations.  In  physical  terms, 
this  means  that  the  response  of  a fluid  parcel  to  a vertical  perturbation  is  associated  with 
changes  in  the  background  density  gradient  rather  with  any  further  compression  of  the 
parcel.  That  is,  although  the  basic-state  density  profile  develops  as  a result  of  the  large 
compressibility  of  the  fluid  near  its  critical  point,  acoustic  waves  do  not  interact  strongly 
with  internal  gravity  waves.  This  can  be  put  in  more  quantitative  terms  if  we  consider  the 
Brunt- Vaisala  frequency 


p dz 


(45) 


where  c is  the  acoustic  sound  speed.  This  quantity  measures  the  fluid’s  oscillatory  response 
to  stratification  and  compression.^'^  Berg  et  al.  argued  that  for  the  near-critical  xenon  system 
under  consideration,  the  Brunt-Vaisala  frequency  could  be  approximated  by  the  stratification 
term  alone.  This  was  based  on  a direct  comparison  for  xenon  of  the  first  and  second  terms 
in  equation  (45).  Since  we  expect  that  our  density  profiles  approximate  theirs  in  the  bulk 
regions  and  that  in  the  interfacial  layer  the  density  gradients  used  here  may  be  quite  large, 
we  anticipate  that  this  is  a reasonable  approximation  to  pursue  here.  This  can  be  stated 
more  formally  in  terms  of  incompressible  perturbations  if  we  consider  the  equation  of  state 
p = p[s,p).  It  follows  that 


Dp 

'm 


dp 


Ds  dp 


Dp  dp  Ds  2^P 
Dt  ds  Dt  Dt 

p 


(46) 


Since  the  base  state  has  no  flow  and  varies  only  in  the  vertical  direction  we  can  reduce  this 
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to  a statement  about  the  perturbed  quantities  (denoted  by  primes) 


dp  , 'dpQ 


' dp  / dpo 


dp 

ds 


ds  idso 


(47) 


In  terms  of  dimensionless  quantities,  the  condition  LgjcP'  <C  1 reveals  at  leading  order  the 
condition  of  incompressible  perturbations  dp  /dt  + w dpo/dz  = 0,  or  V • u =0  using  (la). 

We  insert  expansions  (44)  into  the  governing  equations  (1),  use  the  condition  of  incom- 
pressible perturbations,  linearize,  and  find  that  the  perturbation  quantities  satisfy 


• - , dpo 

lujp  -f  w— — == 

dz 

dw 

q^u  + qyV  — = 

dz 


icopou  - q^p  = 


iiopov  — qyp  = 


dp 

lujpow  4-  — = 
dz 


0, 

0, 


Kpoqx 


d^p\ 

dzy  ’ 


Kpoqy  9 ^ - VT 


-pg  -t-  Kpo 


dz^J 


+ Kp 


d^pQ 
dz^  ’ 


(48a) 

(48b) 

(48c) 

(48d) 

(48e) 


where  q'^  = ql  -h  q^-  Note  that  by  virtue  of  the  incompressible  perturbation  assumption, 
the  equation  for  entropy  s decouples  and  is  not  needed  to  compute  the  frequency.  These 
equations  can  be  comoined  in  a straightforward  manner  to  obtain  the  single  perturbation 
equation 


1 - 

d^w 

■jV2 

^ q^  1 d{poM^) 

dw  , 

— q 

1 - ^ 

[ u;2  J 

dz^ 

. 9 

pq  dz 

dz 

where 


(49) 


9 dpo 
Po  dz  ’ 

K 

Po  \dz  J 


(50a) 

(50b) 


The  boundary  conditions  are  that  the  vertical  component  of  the  velocity  vanishes  on  the 
upper  and  lower  boundaries  w{L)  = w[  — L)  = 0.  This  is  an  eigenvalue  problem  where  the 
eigenvalue  is  the  frequency  io  and  the  eigenfunction  w characterizes  the  wave  mode. 

It  is  of  particular  interest  to  note  that  this  is  a second-order  system.  The  classical  model 
us^'-d  by  Berg  et  al.  for  the  single-phase  region,  which  is  also  second-order,  can  be  obtained  by 
taxing  K = 0 (i.e.  = 0)  in  equation  (49).  That  is,  the  inclusion  of  the  nonclassical  terms 

in  the  incompressible  limit  does  not  result  in  a higher-order  differential  equation  relative  to 
that  obtained  for  the  classical  model.  However,  we  note  here  (and  show  later)  that  if  we 
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include  the  effects  of  compressible  perturbations  the  result  is  a set  of  two  coupled  second- 
order  equations.  We  present  those  equations  and  discuss  them  in  more  detail  in  section  7. 

We  have  solved  the  eigenvalue  problem  (49)  using  a method  given  by  Keller^^  whereby  we 
replace  one  of  the  homogeneous  boundary  conditions  with  an  independent  inhomogeneous 
one.  The  integration  of  the  resulting  inhomogeneous  boundary  value  problem  was  done 
using  SUPORT^®.  This  code  uses  a superposition  of  linearly-independent  solutions  coupled 
with  an  orthonormalization  procedure  to  maintain  their  (numerical)  independence.  The 
integrations  were  performed  with  either  an  Adams-type  method  or  a Runge-Kutta  method. 
The  eigenvalue  u)  and  eigenfunction  w[z)  which  satisfy  the  original  homogeneous  boundary 
condition  were  then  obtained  iteratively. 

Figure  4 shows  the  internal  wave  frequencies  obtained  experimentally  (solid  points)  and 
theoretically  using  the  restricted  cubic  model  (dashed  curves)  by  Berg  et  al.  and  the  theoreti- 
cal predictions  of  the  present  diffuse-interface  approach  using  the  van  der  Waals  model  (solid 
curves).  The  parameter  values  used  for  the  diffuse-interface  calculations  are  cdl  = 10“'^, 
a = 4.85  and  g = 1.631  X 10  The  value  of  e^L  used  here  was  chosen  somewhat  arbitrarily 
but  is  well  within  the  convergence  region  (with  respect  to  the  internal  wave  frequencies). 

If  we  first  focus  on  the  experimental  data,  we  see  that  there  are  two  separate  modes, 
whose  frequencies  have  different  behavior  with  respect  to  the  temperature.  Based  on  the 
eigenfunctions  obtained  theoretically  by  Berg  et  ai,  corresponding  to  the  eigenvalues  shown 
by  the  dashed  curves,  we  know  that  the  upper  mode  corresponds  to  a wave  mode  in  which 
there  is  relatively  large  amplitude  disturbances  near  the  horizontal  mid-plane  of  the  cell. 
Above  the  critical  temperature  this  corresponds  to  an  internal  wave  set  up  over  the  whole 
cell,  while  below  the  critical  temperature,  this  mode  corresponds  to  a sloshing-type  mode 
with  relatively  large  interface  deflection.  The  lower  mode  has  relatively  little  motion  near 
the  center  of  the  cell.  Above  the  critical  temperature,  this  corresponds  to  internal  waves  set 
up  in  both  the  upper  and  lower  portions  of  the  cell,  while  below  the  critical  temperature 
this  corresponds  to  a relatively  quiescent  lower  phase  and  internal  wave  motion  in  the  upper 
phase  only. 

There  are  a number  of  comparisons  that  can  be  made  between  the  theoretical  results 
predicted  with  the  diffuse-interface  model  and  those  obtained  by  Berg  et  al.  using  a classical 
approach.  While  by  design,  the  diffuse-interface  model  agrees  with  the  classical  theory  for 
temperatures  well  above  the  critical  temperature  (recall  that  the  parameters  in  the  present 
equation  of  state  were  chosen  so  that  the  density  profile  matches  that  of  the  Berg  et  al. 
theory  well  above  the  critical  temperature),  the  difference  becomes  more  pronounced  as 
the  temperature  is  decreased.  As  we  shall  see  below  and  in  the  next  figure,  this  difference 
can  be  attributed  in  full  to  the  difference  between  the  two  equations  of  state  used  in  the 
two  models.  That  is,  the  diffuse-interface  results  differ  from  the  classical  ones  of  Berg  et 
al.  because  the  variation  with  temperature  of  the  density  profiles  obtained  with  the  van 
der  Waals  equation  of  state  used  here  does  not  precisely  match  that  of  the  density  profiles 
obtained  using  the  restricted  cubic  model.  Another  result  of  the  diffuse-interface  model  is 
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the  premature  termination  of  the  frequency  curves  below  the  critical  temperature.  As  we 
shall  see  below,  this  is  the  result  of  a breakdown  of  the  disturbance  equation  (49)  as  the 
coefficient  of  the  term  d^wjdz^  vanishes. 

Figure  5 shows  a comparison  of  the  frequency  computed  using  the  diffuse-interface  model 
and  that  computed  using  the  classical  model  of  Berg  et  ai,  where  in  each  case  the  equation 
of  state  (40)  is  used.  Again  the  solid  curves  show  the  diffuse-interface  results  (with  the  same 
parameter  values  as  shown  in  figure  4)  and  the  dashed  curves  show  the  classical  results.  This 
comparison  shows  that  with  the  exception  of  the  regions  below  the  critical  temperature  where 
the  disturbance  equation  breaks  down  (see  below),  the  diffuse-interface  results  reproduce  the 
classical  results.  This  indicates  that  it  is  the  use  of  the  simpler  equation  of  state  (40)  in  the 
present  model,  rather  than  an  inherent  difference  between  the  diffuse  and  classical  models, 
which  accounts  for  the  differences  between  the  theoretical  predictions  shown  in  figure  4. 

The  disturbance  equation  (49)  is  a second-order  equation  whose  leading  coefficient  1 — 
as  noted,  may  vanish.  Although  this  coefficient  depends  on  the  vertical  coordinate, 
we  can  estimate  the  point  at  which  it  first  vanishes  by  evaluating  it  at  i = 0,  where  the 
density  gradient,  and  hence  M,  is  largest.  In  terms  of  the  frequency  = a;/27r  (in  Hz)  and 
temperature  T this  boundary  is  given  by 

n _ {qL)a^  T-T^ 
yJJJl  27rV^  Tc 

Note  that  this  result  does  not  depend  on  the  thickness  of  the  interface.  We  have  plotted 
this  boundary  in  figure  5 (dashed-dotted  line)  for  the  parameter  values  as  given  above. 
This  boundary  is  consistent  with  the  points  at  which  we  can  no  longer  compute  numerical 
solutions  to  the  present  model  (i.e.  where  the  solid  curves  terminate  in  figures  4 and  5). 

7.  Discussion 

In  this  section  we  discuss  several  of  the  issues  raised  by  the  above  described  analysis  and 
numerical  results.  In  particular  we  shall  address  possible  improvements  in  terms  of  the 
equation  of  state  used  in  the  diffuse-interface  model  and  also  the  issue  of  the  limitations  of 
the  incomnressible  perturbation  assumption. 

The  aa vantage  of  the  equation  of  state  employed  in  the  present  analysis  is  that  from  both 
a physical  and  mathematical  viewpoint  it  is  the  simplest  characterization  of  the  near-critical 
behavior.  The  disadvantage,  as  we  have  seen,  is  that  it  does  not  provide  as  accurate  of  a 
representation  of  the  density  profiles  as  does  the  restricted  cubic  model,  for  example.  One 
possible  improvement  which  can  be  made  is  to  use  a modified  van  der  Waals  equation  of  state 
(e.g.  Rowlinson  and  Widom^,  Fisk  and  Widom^^).  The  modified  van  der  Waals  equation  of 
state  improves  on  the  predictions  in  terms  of  critical  exponents  relative  to  the  van  der  Waals 
equation  of  state  but  still  allows  an  analytical  solution  for  the  density  profiles. 


(51) 
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The  numerical  solution  and  analysis  of  the  disturbance  equation  (49)  showed  that  its 
leading  coefficient  can  vanish.  This  indicates  the  presence  of  a singularity.  A key  assumption 
used  to  derive  this  equation  was  that  the  perturbations  were  incompressible.  Although  we 
have  used  a simple  equation  of  state  and  have  used  approximate  techniques  to  write  down 
the  basic-state  density  profile,  it  seems  quite  clear  that  a more  accurate  representation  of  the 
density  profile  shall  ultimately  suffer  the  same  consequences.  In  order  to  further  understand 
this  singularity,  we  shall  rederive  the  disturbance  equations  without  making  the  assumption 
of  incompressibility.  Using  the  same  approach  as  described  in  section  6,  we  find  that  the 
disturbance  equations  that  retain  the  effects  of  compressibility  are  given  by 


. .dp  ^ (Pp  ^ ^ dq  . d?q 

A\p  + A2—  h + A4q  -f  As-  f-  Aq-—  = 0, 

CJ/Z  CLZ  CbZ  CLZ 

Bip  + ^2^  + B^q  Bs^  A Be^  = 0, 
dz  dz  dz^ 


(52a) 

(52b) 


where  q = pow  and  the  (variable)  coefficients  Ai  and  Bi  are  given  in  appendix  C. 

These  are  two  coupled  second-order  equations  describing  a ‘diffuse/compressible’  flow. 
It  is  of  interest  to  recover  several  lower-order  cases  from  this  system;  that  is,  the  ‘dif- 
fuse/incompressible’, ‘classical/compressible’  and  ‘classical/incompressible’  cases. 

The  ‘diffuse/incompressible’  case  can  be  obtained  formally  by  setting  ljp{z)  ==  0 in  the 
expressions  in  appendix  C.  In  this  case,  the  coefficients  A2,  A3  and  B2  all  vanish.  It  can  then 
be  easily  seen  that  equations  (52)  reduce  to  a second-order  equation  for  q,  or  equivalently, 
for  w (see  equation  (49)).  Therefore,  in  the  context  of  the  diffuse-interface  formulation,  the 
incompressible  limit  is  singular. 

This  is  in  contrast  to  the  classical  case,  as  can  be  seen  if  we  compare  the  ‘classi- 
cal/compressible’ and  ‘classical/incompressible’  cases.  First,  to  obtain  the  classical  (single- 
phase) description  from  equations  (52)  we  take  K = 0.  Here,  the  coefficients  A2,  A3  and  Ae 
vanish.  It  is  again  easy  to  see  that  the  resulting  coupled  equations  for  p and  q can  be  reduced 
to  a second-order  equation  for  q.  This  represents  the  ‘classical/compressible’  case.  In  the 
‘classical/incompressible’  case,  which  again  we  can  obtain  by  further  setting  l/P[z)  = 0,  we 
find  that  the  coefficient  B2  vanishes.  Therefore,  in  the  classical  case,  the  disturbance  equa- 
tions are  second-order  regardless  of  whether  or  not  the  effects  of  compressibility  are  included. 
The  above  description  shows  that  the  limit  of  incompressible  perturbations  is  a regular  limit 
in  the  classical  formulation  but  is  a singular  limit  in  the  diffuse-interface  formulation. 

In  order  to  solve  equations  (52),  the  sound  speed  must  first  be  specified.  This,  however, 
cannot  be  calculated  from  the  classical  van  der  Waals  equation  of  state  alone.  A ther- 
mal equation  of  state,  consistent  with  the  van  der  Waals  equation,  must  be  given  (e.g.  see 
Callen33). 

For  the  present  purposes,  we  would  like  only  to  confirm  that  the  inclusion  of  the  terms 
due  to  compressibility  relieve  the  singularity  observed  in  the  incompressible  calculations  and 
allow  us  to  calculate  frequencies  at  lower  temperatures.  Therefore,  we  shall  not  be  concerned 
with  the  specific  form  used  for  the  sound  speed.  In  the  single-phase  region,  we  have  used 
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the  ‘linear  model’  as  given  by  Hohenberg  and  Barmatz^®  to  compute  the  sound  speed.  Note 
that  this  model  is  used  to  compute  the  sound  speed  only;  we  still  use  the  van  der  Waals 
equation  of  state  to  compute  the  density  profiles.  Our  calculations  show  that  the  inclusion 
of  compressibility  effects  in  the  one-phase  region  lead  to  negligible  corrections  in  terms  of 
the  resulting  internal  wave  frequencies.  This  is  to  be  expected  based  on  the  fact  that  our 
incompressible  calculations  encountered  no  difficulty  in  the  one-phase  region.  That  is,  the 
inclusion  of  compressibility  in  the  one-phase  region  represents  a regular  perturbation.  This 
reinforces  the  statements  given  by  Berg  et  al.^^  who  argued  that  the  compressibility  effects 
should  be  negligible  with  respect  to  the  internal  gravity  wave  frequency  predictions.  The 
above  described  ‘linear  model’  applies  in  the  classical  sense  below  the  critical  point,  in  that  it 
can  be  used  to  calculate  the  sound  speed  on  either  side  of  a sharp  interface.  However,  there 
is  no  rigorous  way  of  connecting  the  two  profiles  through  the  diffuse-interface  used  in  the 
present  calculations.  Therefore,  in  the  two-phase  region  we  have  used  a constant  value  for 
the  sound  speed  in  order  to  solve  the  eigenvalue  problem.  Although  this  ad  hoc  prescription 
for  the  sound  speed  does  not  provide  an  accurate  approximation  to  the  true  sound  speed, 
which  can  vary  significantly  in  the  vertical  direction,  it  is  sufficient  for  our  purposes  in  that 
it  is  nonzero  everywhere  and  can  be  made  to  represent  the  true  sound  speed  in  an  average 
sense.  We  have  found  that  this  approach  does  in  fact  allow  us  to  calculate  internal  wave 
frequencies  beyond  the  line  of  singularity  shown  in  figure  5.  Therefore,  we  can  conclude  that 
the  effects  of  compressibility  relieve  this  singularity  and  may  not  be  negligible  when  using 
the  diffuse-interface  formulation  as  we  have  presented  here. 

As  an  additional  check,  we  have  calculated  the  solutions  for  the  case  where  the  capillarity 
terms  are  included  in  the  basic-state  density  profile  while  the  terms  involving  and  its 
derivatives  are  ignored  in  an  ad  hoc  way  in  equation  (49).  From  the  point  of  view  of  the 
perturbation  equation,  the  system  behaves  as  if  it  were  in  a single  phase,  with  the  density 
profile  represented  as  before  in  both  the  one-phase  and  two-phase  regions.  The  results  of  this 
calculation  are  graphically  identical  to  the  dashed  curves  shown  in  figure  5,  which  represent 
the  sharp-interface  model  results.  Unlike  the  full  incompressible  case  where  the  capillary 
terms  are  included  in  the  perturbation  equation  (49),  these  calculations  are  not  limited 
by  the  line  of  singularity. 

8.  Conclusions 

We  have  presented  diffuse-interface  equations  for  an  inviscid,  compressible  and  adiabatic 
flow  and  have  described  these  equations  in  terms  of  global  quantities  and  balance  laws.  The 
key  difference  between  the  diffuse-interface  formulation  and  the  classical  approach  is  the 
presence  of  a capillary  tensor  in  the  momentum  balance.  This  tensor  models  the  capillary 
forces  associated  with  the  diffuse  interface.  In  contrast  to  the  classical  formulation,  no 
interface  needs  to  be  tracked. 

We  have  analyzed  these  equations  and  found  that  we  can  recover  the  classical  equations 
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and  interfacial  jump  conditions  for  a fluid-fluid  system  in  the  sharp-interface  limit.  The 
interfacial  tension  is  found  to  be  consistent  with  the  excess  Kramer’s  (Grand  canonical) 
potential  energy. 

We  have  used  this  diffuse-interface  model  to  describe  internal  gravity  waves  in  near- 
critical  xenon.  This  analysis  is  an  extension  of  that  done  by  Berg  et  al.^^ , who  studied  the 
problem  experimentally  and  theoretically  using  classical  hydrodynamic  descriptions  above 
and  below  the  critical  temperature.  The  present  diffuse-interface  description  of  the  hydro- 
dynamics can  be  used  to  calculate  the  internal  wave  frequencies  both  above  and  below  the 
critical  temperature.  Our  results  show  that  the  diffuse-interface  model  recovers  the  classical 
results  when  the  same  equation  of  state,  used  to  compute  the  equilibrium  density  profiles,  is 
used  in  each  case. 

An  unexpected  result  identified  here  is  that  the  incompressible  flow  limit  is  singular  in 
the  context  of  the  diffuse-interface  limit.  That  is,  when  the  effects  of  compressibility  are 
neglected,  the  internal  wave  problem  is  governed  by  a second-order  perturbation  equation, 
but  when  the  effects  of  compressibility  are  included,  the  problem  is  governed  by  a system  of 
two  second-order  perturbation  equations.  This  is  in  contrast  to  the  classical  case  where  the 
limit  of  incompressibility  is  regular;  both  cases,  compressible  and  incompressible,  result  in  a 
second-order  perturbation  equation  describing  the  wave  modes. 
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Appendix  A.  Integral  Quantities 


The  following  expressions  are  used  in  calculating  the  time  rates  of  change  of  the  total  mass, 
momentum,  energy  and  entropy.  First,  we  note  that  for  a material  volume  which  moves 
with  the  fluid  the  identity 


dt 


[ 4>dV=[  (^  + y-i^u)]dV, 

jQ{t)  Jn{t)  \dt  J 


(Al) 


holds  for  the  scalar  quantity  (f). 

Note  that  the  following  results  make  use  of  the  condition  of  mass  conservation  dM/dt  — 0, 
or 
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Appendix  B.  Dissipative  Effects 

Here  we  formulate  the  diffuse-interface  equations  that  account  for  thermal  and  viscous  dissi- 
pation. Our  approach  follows  the  formalism  used  in  notes  by  Sekerka^®  on  entropy  produc- 
tion. 

We  start  with  the  same  total  quantities  for  mass  M and  energy  E as  used  in  Section  2 with 
the  exception  that  the  gradient  energy  coefficient  is  now  assumed  to  depend  on  temperature 
T . The  total  entropy  S now  includes  both  classical  and  nonclassical  contributions  and  a 
temperature-dependent  gradient  entropy  coefficient.  Additionally,  we  include  a definition 
for  total  momentum  P . These  quantities,  defined  on  the  material  volume  0(i),  are 

M = f pdV,  (Bla) 

Jn{t) 
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p = 

f pudV, 

Jn{t) 

(Bib) 

E = 

(Blc) 

S = 

(Bid) 

First,  we  require  that  the  mass  in  any  subvolume  Q,{t)  is  conserved,  giving  dM j dt  = 0 

and 

+ = 0.  (B2) 

Next,  we  write  down  the  general  forms  of  physical  balance  laws  associated  with  P,  E 
and  5, 
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where  m is  a general  stress  tensor  which  may  include  both  classical  and  nonclassical  stresses, 
q is  the  classical  conductive  heat  flux,  is  a nonclassical  energy  flux,  is  a nonclassical 
entropy  flux  and  is  the  volumetric  rate  of  entropy  production  which  may  include  both 
classical  and  nonclassical  contributions.  In  physical  terms,  the  momentum  balance  (B3a) 
states  that  changes  in  momentum  result  from  forces  on  the  boundary  and  body  forces  in 
the  interior  (here  we  are  considering  gravity  as  the  only  such  body  force).  The  energy 
balance  (B3b)  states  that  changes  in  energy  result  from  the  rate  of  work  done  by  the  forces 
on  the  boundary,  the  classical  heat  flux  and  a nonclassical  energy  flux  through  the  boundary. 
The  entropy  balance  (B3c)  states  that  the  rate  of  change  of  entropy  in  the  control  volume 
minus  the  fluxes  of  entropy  (classical  and  nonclassical)  through  the  boundary  must  be  equal 
to  the  entropy  production.  The  Second  Law  of  Thermodynamics  requires  that  this  entropy 
production  be  positive  for  dissipative  systems. 

In  what  follows,  we  use  the  definitions  (Bl)  to  derive  identities  for  dP jdt,  dE/dt  and 
dSjdt  with  the  aid  of  the  results  given  in  Appendix  A.  We  then  equate  these  expressions  with 
the  balance  laws  (B3)  to  obtain  local  balance  laws  in  terms  of  the  general  quantities  m,  g, 
3^e^\  Finally,  we  appeal  to  the  Second  Law  of  Thermodynamics  > 0) 

to  identify  thermodynamically  consistent  forms  for  m,  g,  and 

Momentum:  We  seek  a local  balance  law  for  momentum  by  calculating  the  time  rate 
of  change  of  the  total  momentum  given  by  equation  (Bib)  and  then  equating  the  result  with 
the  balance  law  for  momentum  (B3a).  We  use  the  divergence  theorem  to  express  boundary 
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integrals  as  volume  integrals  and  note  that  the  domain  of  integration  Q(t)  is  arbitrary.  This 
gives  the  local  momentum  balance 


Du 


V • m + pgz  = 0. 


(B4) 


Energy:  Next,  in  the  same  manner  as  described  above,  we  obtain  the  local  energy 
balance 


/ Du  „ De  „ „ „ 
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Here  we  have  defined  T = (/9V^/?+^|V/jp)J— V/O0  Vp,  which  we  can  identify  as  a Korteweg- 
like  capillary  tensor^'^.  We  shall  find  below  that  the  tensor  m is  related  to  T.  We  can  simplify 
the  local  energy  balance  by  subtracting  the  contribution  due  to  mechanical  energy,  which 
is  expressed  as  the  velocity  u dotted  with  the  momentum  balance  (B4).  This  leads  to  the 
simplified  local  energy  balance 
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We  can  express  this  equation  in  terms  of  the  entropy  density  s by  first  noting  that  the  energy 
density  e is  a function  of  s and  p and  then  appealing  to  the  thermodynamic  relationships 
dejds  — T and  de/dp  = pj p^  which  allow  us  to  write 


De  ^ rj,Ds  p Dp 
Dt  Dt  p'^  Dt 


Combining  this  with  equations  (B2)  and  (B6)  yields 


pT 


Ds  Dp 

+ V • g = pV  u m : Vtt  — Ke{T)T  ; Vit  + • Vp 

--\Wp\^P^  -V  ■( 


- V - {Ke(T)PAVp  + 


(B7) 


(B8) 


Entropy:  Finally,  we  can  obtain  a local  entropy  balance  using  equations  (Bid)  and 
(B3c).  This  equation  is  given  by 

4 + ^ = -Ks{T)T:Vu  + ^VKs-Vp-\\Vp\^P§^ 

-q  ■ V (^)  + - V . {Ks(T)^Vp  + • (B9) 

We  assume  that  m can  be  written  as  the  sum  of  a purely  classical  part,  which  we  take  to 
be  the  stress  tensor  cr  = —pi  + t where  t is  the  deviatoric  stress  tensor,  and  a nonclassical 
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part  so  that  m = <t  + As  a consequence  of  this  and  equations  (B8)  and  (B9) 

we  obtain  the  following  expression  for  the  volumetric  rate  of  entropy  production 


^prod 


aV|- 


r : Vu  ^ 1 

^ T ^ f 


- [Ke{T)  - TKs{T)]t)  ; Vu 


+ ^ (VA:^  - TVifs)  - 

Dt  ^ ^ 2'  T \ Dt  Dt  J 

+V  ■ + - iv . (ifE(r)|^Vp+iM) . (bio) 


We  can^^simplify  the  expression  (BIO)  using  the  following  argument  which  relates  the 
gradient  energy  and  gradient  entropy  coefficients.  Suppose  we  write  the  total  free  energy 
density  (and  similarly  e‘°‘  and  5*°‘)  as  a function  of  not  only  of  the  classical  variables 
but  also  a = as  defined  below 


f^\T,p,a)  = 

f{T,p)  + KF{T)^-\Vp\\ 

(Blla) 

e°\s,p,a)  = 

e{3,p)  + KE{T)^\Vp\\ 

(Bllb) 

s^°\T,p,a)  = 

s(T,p)  + Ks(T)^-\Vp\\ 

(Bile) 

Further,  we  assume  that  the  same  thermodynamic  relations  which  hold  for  the  bulk  quantities 
also  hold  for  the  total  quantities.  In  particular,  we  assume  = e*°‘  — Ts‘°*  and  = 
—fj^.  Here  the  subscript  denotes  differentiation  with  respect  to  temperature  while  the  other 
variables  (i.e.  p and  a)  are  held  fixed.  Then,  if  we  use  f = e — Ts,  s = — /t  we  must  also  have 
the  relations  Kf{T)  = Ke{T)  - TKs{T),  Ks[T)  = -K'e[T),  Ke{T)  = Ke{T)  - TK'e{T)  so 
that  K'g{T)  = —Ke{T)  and  K'^[T)  = —TKe{T)  = TKg{T)  where  ' denotes  the  derivative 
with  respect  to  temperature.  These  results  lead  to  a simplified  expression  for  the  volumetric 
rate  of  entropy  production 


= g . V (i)  + 1 (m(“)  - : V-u 

+V  . {ks{T)^Vp  + it'')  - iv  . (^Ke{T)^Vp  + jt^)  . (B12) 

Our  next  objective  is  to  identify  forms  for  the  quantities  q,  r,  and  which 

are  consistent  with  positive  entropy  production. 

Since  equation  (B12)  must  hold  when  the  nonclassical  terms  are  absent  we  can  first  focus 
on  the  classical  terms.  In  the  classical  case  we  can  guarantee  positive  entropy  production 
locally  by  choosing 


q = kT^V 


(B13) 


and 


r = p,{Vu  + Vtfc^) /i(V  • u)J, 


(B14) 
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where  k is  the  thermal  conductivity  and  //  is  the  dynamic  viscosity  of  the  fluid.  The  form 
for  q is  the  Fourier  Law  for  heat  conduction  (e.g.  see  Kittel  and  Kroemer'^°  or  Carslaw  and 
Jaeger"^^).  The  form  for  the  deviatoric  stress  tensor  r is  just  that  for  a Newtonian  fluid 
(e.g.  see  Batchelor'^^). 

We  now  consider  the  nonclassical  contributions  to  the  entropy  production  (B12).  This 
form  suggests  that  the  surface  tension  is  associated  with  the  excess  free  energy,  which  is 
consistent  with  what  we  expect  based  on  thermodynamic  arguments  (see  Section  3).  Further, 
we  make  the  ansatz  that  the  effects  of  surface  tension  are  reversible  in  nature.  This  implies 
that  if  the  effects  of  classical  viscous  and  thermal  dissipation  are  neglected,  the  entropy 
production  is  zero.  The  following  specifications  for  and  j^s^\ 


(nc) 

‘ = 

Kf{T)T, 

(B15a) 

3e  — 

-KdT)fvp, 

(B15b) 

3 s — 

-Ks(T)^Wp, 

(B15c) 

are  consistent  with  these  conditions.  This  says  that  the  nonclassical  stress  tensor 
corresponds  to  the  Korteweg  capillary  tensor  T with  the  coefficient  associated  with  the 
gradient  free  energy.  The  nonclassical  flux  terms  are  in  the  direction  of  the  density  gradient 
and  are  zero  if  the  flow  is  incompressible  (i.e.  Dp/Dt  = 0).  These  nonclassical  energy 
and  entropy  fluxes  are  associated  with  the  gradient  internal  energy  and  gradient  entropy, 
respectively.  A similar  nonclassical  entropy  flux  term  was  identified  in  the  phase-field  model 
of  solidification  derived  by  Wang  et  al?’^  (see  their  equation  (6)).  Whereas  their  term  involved 
the  partial  derivative  of  the  order  parameter  with  respect  to  time,  our  term  involves  the  total 
derivative,  since  we  have  accounted  for  fluid  motion.  They  identified  this  term  as  an  entropy 
flux  associated  with  variations  in  the  phase-field  at  the  boundary  of  the  subvolume. 

The  above  prescription  leads  to  the  set  of  governing  equations 

= — pV  • It,  (B16a) 

= V - a - pgz  + V ■{Kf{T)T),  (B16b) 

= CT-.WU-Wq-  TKs(T)T  : Vu  + ^VKe  ^ Vp  - 
= r : V«  - V ■ g - TKs(T)T  : V«  + ■ Vp  - ilV^|";^.(B16d) 

V ^ (Kf(T)T)  = pV  [Kf{T)V^p)  - ]^\Vp\'‘VKe(T).  (B17) 


Dt 

Du 

De 

rr^Ds 


Note  that 
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Appendix  C.  Coefficients  of  Compressible  Equations 


The  following  are  the  coefficients  appearing  in  the  compressible  perturbation  equations  (52). 
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2D  view: 


Figure  1:  The  pillbox  enclosing  the  interfacial  region  shown  in  a side  view. 
We  identify  the  constant  density  surface  that  is  the  limiting  surface  of  the 
interfacial  region  as  its  thickness  goes  to  zero.  This  surface  has  the  normal 
vector  n.  The  pillbox  has  normal  vectors  iitop,  Ubot  and  m on  its  top,  bottom 
and  side,  respectively.  The  figure  also  illustrates  the  limiting  case  considered 
in  the  pillbox  argument.  The  upper  and  lower  surfaces  of  the  pillbox,  defined 
by  C = where  is  a local  normal  coordinate,  are  squeezed  together  in  such  a 
way  that  the  variation  in  density  p through  the  interfacial  region  is  completely 
enclosed  within  the  pillbox.  That  is,  we  consider  the  limit  e <C  <C  T where 
e is  a measure  of  the  interfacial  thickness  and  L is  an  outer  0(1)  length  scale. 
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Figure  2:  This  figure  shows  the  model  configuration  filled  with  a stratified 
fluid.  The  dimensions  are  0 < x < Ui,  0 < y < Uy  and  —L  < z < L where 
Ui  = 7.6mm,  ay  = 38mm  and  L = 9.5mm  are  the  values  corresponding  to  the 
experiment  by  Berg  et  al.  The  paddle,  which  in  the  experiments  creates  the 
initial  disturbance,  is  shown  for  reference  and  is  not  present  in  the  mathemat- 
ical model. 
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Figure  3:  This  figure  shows  the  static  density  profiles  representing  solutions  of 
equation  (41).  Each  curve  corresponds  to  a different  temperature  as  indicated. 
The  dashed  curve  corresponds  to  the  density  profile  at  the  critical  tempera- 
ture. Above  the  critical  temperature,  there  is  a single  stratified  phase.  Below 
the  critical  temperature,  the  fluid  separates  into  two  stratified  phases.  The 
parameter  values  used  to  calculate  these  profiles  are  e^L  — 10“^,  a — 4.85, 
g = 1.631  X 10-A 
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T-Jc  (mK) 


Figure  4:  This  figure  shows  the  internal  wave  frequencies  Q = uj{2'k  (in  Hz)  ob- 
tained experimentally  by  Berg  et  al.  (data  points),  the  theoretical  predictions 
by  Berg  et  al.  (dashed  curves)  using  two  separate  models  above  and  below  the 
critical  temperature  coupled  with  a restricted  cubic  model  for  the  equation  of 
state,  and  the  theoretical  predictions  of  the  present  diffuse-interface  approach 
(solid  curves)  using  a van  der  Waals  equation  of  state.  The  vertical  dashed 
line  indicates  the  critical  temperature.  The  parameter  values  used  for  the 
diffuse-interface  calculations  are  em  = a = 4.85  and  g = 1.631  x 10“'^. 
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T-T^  (mK) 


Figure  5:  This  figure  shows  a direct  comparison  of  the  frequency  calculated 
using  the  classical  model  and  that  calculated  with  the  diffuse-interface  model. 
In  each  case  we  have  used  (40)  as  the  equation  of  state.  The  vertical  dashed 
line  indicates  the  critical  temperature.  We  find  that  in  this  case  the  diffuse- 
interface  model  (solid  curves)  recovers  the  classical  model  (dashed  curves) 
given  by  Berg  et  al.  The  dashed-dotted  line  in  this  figure  indicates  the  location 
where  the  coefficient  1 — q^M'^  fuj'^  of  the  second-order  term  in  equation  (49)  is 
predicted  to  vanish,  and  hence  indicates  the  boundary  beyond  which  we  can- 
not compute  with  the  present  model  in  the  incompressible  perturbation  limit. 
Since  this  coefficient  depends  on  the  vertical  coordinate,  we  have  identified  this 
boundary  by  evaluating  the  coefficient  at  i = 0,  where  the  density  gradient, 
and  hence  M,  is  greatest. 
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